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ABSTRACT 

We present an equilibrium statistical mechanical theory of collisionless self-gravitational systems 
with isotropic velocity distributions. Compared to existing standard theories, we introduce two 
changes: (1) the number of possible microstates is computed in energy (orbit) space rather than phase 
space and (2) low occupation numbers are treated more appropriately than using Stirling's approxi- 
mation. Combined, the two modifications predict that the relaxed parts of collisionless self-gravitating 
systems, such as dark-matter halos, have a differential energy distribution N{e) oc [exp(</>o — e) — 1], 
dubbed "DARKexp". Such systems have central power-law density cusps p(r) oc r _1 , which suggests 
a statistical mechanical origin of cusps in simulated dark-matter halos. 

Subject headings: dark matter — galaxies: halos — galaxies: kinematics and dynamics — methods: 
analytical 



1. INTRODUCTION 

The apparent universal light distributions in ellipti- 
cal galaxies with two-body collision relaxa tion time ex- 
ceedin g the age of the un i verse motivated lOgorodnikovl 
(|1957f) and ILvnden-Bell (j!967D to seek a fast relax- 
ation process driving the systems toward equilibrium 
in a statistical me chanical sense. In a seminal paper, 
Lvndcn-Bcll ( 1967) introduced the process of 'violent re- 
laxation' (collective energy exchange between rapid po- 
tential fluctuations and individual particles) as respon- 
sible for a short time scale for relaxation. In the same 
paper, he introduced a new kind of statistical mechan- 
ics, that of distinguishable particles subject to an exclu- 
sion principle because collisionless dynamics precludes 
two parcels of phase space from being superimposed. In 
the non-degenerate limit the theory predicts that isother- 
mal spheres are the maximum e ntropy equilibrium st ates 
of the process, as also found by lOgorodnikovl (|1957fl . 

The predictions of the theory are, however, not entirely 
satisfactory. It predicts infinite mass systems despite 
being derived under the constraints of fixed energy and 
mass. It also predicts mass (phase-space density) segre- 
gation despit e the dynamics being col li sionless (i.e., mass 
indep e ndent: iSeverne fc Luwell [19861: iHiorth fe Madsenl 
19931 iNakamural l2000t lArad fc Lvnden-Belll 120051: 
Arad fc Johansson! I2005D . The resulting isothermal 
sphere profile does not reproduce the light profiles of 
elliptical galaxies which was the original motivation. 
There is an arbitrariness in defining the initial states and 
whether to use a par t icle or pha s e element approach (IShul 
19781: iMadsenl 119871: IShul 119871: ISTorth fc Madsenl 119931: 
Kull et all I1996L 119971: IBindoni fc Seccol 12008ft . which 



makes it difficult to assess whether the degenerate limit 
may be relevant. And finally, it is not obvious how to 

1 Dark Cosmology Centre, Niels Bohr Institute, University of 
Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen 0, 
Denmark; jens@dark-cosmology.dk 

2 Department of Astronomy, University of Minnesota, 
116 Church Street SE, Minneapolis, MN 55455, USA; 
llrw@astro.umn.edu 

3 Institute for Theoretical Physics, University of Zurich, Win- 
tcrthurcrstrassc 190, CH-8057 Zurich, Switzerland 



extend the statistical mechanical approach to (spherical) 
systems with anisotropic velocity distributions. 

Some of these issues have been addressed in 
terms of modifications of the extent of the relax- 
ation process, so-called incomplete violent relaxation 
Stiavelli fc Bertinl ll987fl . relaxation in a finite volume 



,Hiorth & Madsen 1991), or explicit scattering processes 
(|Spergel fc Herncmistl 11992). Another approach has 
been to propose a change in the entropy functional 
to be optimi z ed, applicable to non-ex tensive systems 
(|Tsallisl 119881 : iPlastino fc Plastinol 11993ft. although this 
has been demonstrate d not to work (Barnes et al J [20071 : 
iFeron fc Hiorthll2008f ). 

To addr ess the infinite-mass problem, [ Madsen (1996) 
(and later iMenon et ai1l2006l: iDubev et al.ll2008D showed 
that appropriately dealing with small occu pation num - 
bers leads to finite-mass systems, similar to iKingl ([19661 ) 
models, suitable for the description of globular clusters, 
which are driven by collisional relaxation. 

The subject of the origin of universal collisionless self- 
gravitating structure has gained renewed attention in re- 
cent years with the demonstration that the end products 
of numerical simulations of cosmological structure for- 
mation and dark-matter halos have remarkably universal 
profiles, in density as well as in pseudo phase-space den- 
sity, p/<T 3 . A new aspect, not foreseen by Lynden-Bcll's 
theory or any modifications thereof, is that simulated 
dark-matter halos appear to have central density cusps, 
p{r) ~ r _1 , fal ling to r~ 3 or r~ 4 , or even steeper, in 
the o uter parts (iNayarro et al.lll997l 120041 : iMerritt et all 
[20051 INavarro et alJ^OloTT The o rigin of such cusps is 
unclear (see, e.g.. lHenriksenll2006[ ). There also appears 
to be a relation between the local den sity slope and the 
degre e of radial velocity anisotropy ([Hansen fc Moorel 
[20061) . 

In this paper, we explore a new avenue for implement- 
ing the collisionless nature of dark-matter halos. We 
suggest to partition state space in energy-per-unit mass 
shells rather than phase-space elements. Moreover, we 
implement the corresponding effect of small occupation 
numbers in this approach. We show that the resulting 
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state is a finite-mass system with a cusp, r^ 1 , falling in 
the outer parts to r~ 4 for a completely relaxed, isolated, 
isotropic system. This captures the overall properties 
of simulated halos. In companion papers we compare 
in detail the resulting states with numerically generated 
systems suitable for testing our predictions. 

2. STATISTICAL MECHANICS 

We start out, following IBoltzmannl (|1896f) and 
Ogorodnikov (1957), by defining the number of possible 
states: «; 

W = N\Y\ !i± T (1) 

(as discussed in Sectio n 4 we do not introduce an 
exclusion principle, cf . iLvnden-Beil Il967t iShul 119781 : 
IStiavelli fc Bertinlll987[ ). Here, n-i is the occupation num- 
ber in cell of size gi in state space (/i space) . 

It is customary to take the continuous limit (n^ — > ri), 
identify T(n + 1) = n\, and optimize hiW under con- 
straints of fixed total energy and total number of parti- 
cles using a variational approach. This yields 



In g - ip(n + 1) = a + (3E, 



(2) 



where ip(n) = d In T/dn is the digamma function, a and /? 
are Lagrange multipliers associated with the constraints 
of total number of particles and total energy, respectively, 
and E — v 2 /2 + 3>(x) is the energy per unit mass, where 
<i>(x) is the gravitational potential. 
Next, one introduces Stirling's formula 



Inn! = 



I) lnn-n+iln(27r)+^ ; < 8{n) < 1 
2/ 2 12n 

(3) 
(4) 



and uses Stirling's approximation for large n 

In n\ = n In n — n ; n 1 
which implies 

V>(n + l)=lnn; n » 1. (5) 

This leads to 



ff exp(-a-^£). 



(6) 



This is the classical finding for the non-degenerate limit, 
as we did not introduce an exclusion principle in Equa- 
tion (1). In particular, identifying fi space with the 6- 
dimensional (x,v) phase space, the isothermal sphere is 
retrieved if phase space is divided up equally into equal- 
size phase-space cells g: 



/(x,v)=Aexp -/? 



1 



<D(x) 



(7) 



Below we introduce two modifications to this standard 
approach, which, combined, give dramatically different 
structures, which turn out to be reminiscent of the end- 
products of numerical simulations of collisionless self- 
gravitating A-body systems. 

2.1. Orbit Space Versus Phase Space 

The first major modification consists in noting that 
in an equilibrium collisionless system all particles retain 
their energies. This is not the case in an equilibrium colli- 
sional system, such as a classical gas, or a self-gravitating 



system dominated by two-body encounters, such as a 
globular cluster. The fundamental property of a colli- 
sionless system is that particles are distributed on orbits. 
In a relaxed system, once an energy is assigned to a par- 
ticle it stays in a restricted portion of phase space, an 
energy cell. For an isotropic system, energy is the only 
isolating integral. Hence, we argue that in partitioning 
state space, the fundam ental partition is energy spa ce, 
not phase space (see also lEfthvmiopoulos et al.ll2007l) . 

This implies that the occupation number n, Equa- 
tion ([6]), should be interpreted as an indicator of the 
number of particles with a given energy, i.e., 



N(E) oc exp(-/3E) 



(8) 



and not as the number of particles in a parcel of classical 
p hase spa ce , as is usually assumed. 

19821) ( see als o, e.g.. iHernquistl 

ist| r 



Binncv 
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1990; 



Hernquistl 119921 : iBinnev fc Tremaind 11987 ) 
found that elliptical galaxies obeying the i? 1 / 4 law have 
energy distributions very similar to Equation ([8]), how- 
ever these authors' motivation for using this form was 
largely empirical. 

Given N(E) one c an, of course, recover th e classical 
distribution function. IStiavelli fc Bertinl (1987) point out 
that the occupation numbers in the classical phase-space 
can be related to those in the energy space if the former 
are assigned non-equal a priori weights, inversely propor- 
tional to the volume of phase space with a given energy, 



f(E)^g(E)- 1 cM^PE), 



(9) 



from which Equation ([8]) is retrieved since N(E) = 
f(E)g(E) for isotropic systems. The relation be- 
tween the phase-space density / and the differential en- 
ergy distribution N is obtained assuming isotropy from 
/(x, v)d 3 vd 3 x — N(E)dE and the den sity of st ates 
a(E) = d 3 vd 3 x/dE = 16^ 2 J rmax(B) ^/2{E-<f>)r 2 dr 
(jBinnev fc Tremaind fl987T ). 

Ano ther inte r esting aspect of Equation ((8]) , as pointed 
out by iBinnevI ()1982D . is that for systems with less mass 
at the centers than in the outer parts, /3, and hence the 
effe ctive temperature, T = /S^ 1 , need to be negative (see 
also lMerritt et af1 ll"989). This is consistent with the self- 
gravitating nature of the systems, as these are charac- 
terize d by negative heat capacity (|Binney fc Tremainel 
[1981 . 

2.2. Small occupation numbers 
2.2.1. Integer approach 

The second major modification is necessary because 
in systems with a finite potential depth the occupation 
numbers in energy cells with very bound energies can be- 
come small (Equation ((HJ). A similar problem is encoun- 
tered when the classical phase-space distribution func- 
tion, f(E), takes on an exponential form (Equation ([7])), 
but in this case, because the temperature is positive, the 
small occupation numbers are found at the nea r escape 
energi es, i.e., in the outer regions of systems. iMadsenl 
(1996) pointed out that in this case the Stirling approx- 
imatio n, Equation (|4]), breaks down. Following iSimonsI 
(1994), he argued that the appropriate form for Equa- 
tion © is 

Hi = [gt exp(-a - 0Ei)], (10) 
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Figure 1. Approximations to rp(n + 1). Comparison of the stan- 
dard approach = 0, dashed curve) and the approximation used 
here (f = exp(— 7), solid curve). 



where [•] means rounding down to the nearest integer. 

Because the majority of particles in this latter case 
have energies near escape energies, the small occupa- 
tion number modific ation has a dr amatic effect on the 
resulting structures. [Madsenl (|1996| ) showed that Equa- 
tion (|10|) introduces a c utoff which r esults in finite-mass 
systems, similar to the IKind (11966D models of globular 
clusters. In eff ect. iMadsenj (|1996|) analytically derived 
the well-known IKind (|1966l) models, which were origi- 
nally obtained as a simple heuristic modification of the 
isothermal sphere's distribution function, Equation ((7|)|3 

For collisionless systems, where the temperature is neg- 
ative, the cutoff occurs at energies close to that of the 
central potential value, i.e., close to the center of the 
system. Hence, the effect on the structures is not so dra- 
matic in terms of total mass. But as we show below, it 
determines the inner density profile of the equilibrium 
systems. 

2.2.2. Continuous approach 

Because a physical system is not expected to have a 
step-like differential energy distribution N(E), a contin- 
uous version of Equation (fT0|) is needed. We do this by 
introducing a superior approximation to Stirling's for- 
mula which, unlike Equation ([5]), is not limited to large 
occupation numbers. 

We start by noting that when n = 0, one has the ex- 
act result "0(1) = —7, where 7 w 0.57721566... is Euler's 
constant. Combining this with the large n limit, Equa- 
tion ([5]) leads to the approximation 



ip(n + l) wln(n + C)> 



(11) 



with £ = exp(— 7) w 0.561459.... While simple, this ex- 
pression is a remarkably good approximation as shown 
in Figure 1 along with the classical large n approxima- 
tion (Equation ([5]); ( — 0). The deviation from the exact 
ip(n + l) is always positive, with a single maximum differ- 
ence of 0.0237 at n = 0.680. Using this in Equation ©, 
one obtains a modification to Equation ([6]). 



5 exp(-a- (3E) - (. 



(12) 



4 One can possibly improve upon [Madscn (1996)'s Figure 1 by 
assigning different central potential values of the new models that 
match the old models. Then, the difference between the two will 
be confined to low density regions, where it belongs. 



( Landsberg 


1954; 


Dubcv ct al. 


2008) 



The Lagrange multiplier a is determined by requiring 
that n = at some energy E' . Thus, we get 

n = C(exp(-/3[E-E'])-l), (13) 

thereby eliminating the cell size, g, which does not have a 
unique, physically meaningful value in gravitational sys- 
tems. 

For a system which vanishes at the escape energy, 
E' = 0, and so the classical phase-space occupation num- 
ber function is n = £(exp(— f3E) — 1). The systems we are 
interested in vanish at a finite central potential, E' = $o, 
and so the energy-space occupation number function be- 
comes n = ((exp(—(5(E — $o)) — !)■ This is the contin- 
uous version of Equation (fTU| . 

2.3. DARKexp models 

Equation (|13[) with E' = $o incorporates the two mod- 
ifications we introduced in Sections 12.11 and 12.21 Iden- 
tifying the occupation number n as being proportional 
to N, the differential energy distribution, one obtains 
N(E) = A(exp(-p[E — $o\) ~ 1)' where A is determined 
by the mass of the system. In this expression, E and $o 
have units of energy, and /3 is the inverse of temperature. 
We convert this to dimensionless form, using e = j3E and 
4>o — /3$o (note that both e and 0o are positive quantities 
for bound systems), 



N{e) = A{exp[<f> - e] - 1). 



(14) 



We dub this expression the DARKexp. It represents our 
prediction for fully relaxed, collisionless, self-gravitating, 
isotropic systems. Because the arguments presented in 
this paper do not apply to non-isotropic systems, Equa- 
tion (fT4")l cannot be directly compared to the results of N- 
body simulations. Having said that, we note that N(e) 
is not very sensitive to anisotropy, as it depends primar- 
ily on the density profile and not on the dynamics of the 
system. 

The detailed structure of the DARKexp systems, 
including their density and velocity dispersion pro- 
files, will be considered in an accompanying paper 
(Willi ams &: Hjorthl 120101 ). In the next section, we dis- 
cuss the limiting behavior of DARKexp models at small 
and large radii. 

3. RESULTING STRUCTURES 

3.1. Limiting Power-Law Behavior at Small Radii 

In a general spherically symmetric structure, the lim- 
iting form for the central potential can be assumed to 



be 



00 - <t> a r a + 



(15) 



Possion's equation V 2 = AirGp yields a limiting power- 
law behavior of the central density, 

p(r)(xr a - 2 . (16) 

For e —> (f>Q the distribution function then becomes 

/(e) ex (0o - £ )- (4+ " )/2Q (17) 
and the density of states becomes 

ff (e)cx(0o-£) (6+a)/2a . (18) 
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For simi l ar exp ressio ns. seelHenriksen fc Widrowl (|1995() . 
iWidrowl (120 00 ), and |AradetaL| (|2004j). The differential 
mass distribution N(e) — f(e)g(e) then becomes 

N(e)^(cb -e) 1 ^. (19) 

In general, the value of a ranges from for the singular 
isothermal sphere to a = 1 for Navarro-Frenk- White or 
Hernquist profiles, to shallower slopes, reaching a flat 
core for a = 2. Increasing a corresponds to increasingly 
steeper lniV(e) versus e curves. For a — > oo, the system 
develops a hole in the central density profile, which is 
unphysical. 

3.2. Limiting Power- Law Behavior at Small Radii for 
DARKexp Models 

In the DARKexp model, Equation dTJJ, N(e) oc (<po - 
e) as e — > 4>q. In other words, a — 1 and we retrieve the 
central density cusps 

pMocr" 1 , (20) 

known from numerical simulations. The corresponding 
distribution function is 

/(e) cx (0 O - £ r 5 /2 (21) 

for e — » <fio , which is the same as that of the iHernquistl 
(1990) model. 

Thus our proposed differential energy distribution, 
the DARKexp, naturally predicts that the central den- 
sity slopes of collisionless self-gravitating systems should 
asymptote to —1 at the centers of structures. Note that 
this slope is not the result of any specific dynamical 
process operating during the formation of halos, but a 
generic consequence of full relaxation. 

The Appendix discusses the limiting behavior for more 
general cutoffs to N(e) cx exp(— e). 

3.3. Limiting Behavior at Large Radii 

While we are focusing on the inner parts of halos in this 
paper, we note that a full model for the entire system can 
be obtained by assuming that N(e) is finite at the escape 
energy and zero above (as plotted in Figure 2). The ra- 
tionale behind this is that during violent relaxation the 
escape energy is no special location and N (e) is expected 
to be continu ous at what will eventually become the es- 
cape energy (IJaffel IT987t [Hiorth fc Madsenl [l99l . For 
such a model f(e) oc e 5 / 2 and p{r) cx r~ 4 for e — > and 
r — > oo, similar to the Hernquist model and broadly con- 
sistent with numerical simulations of dark-matter halos. 
However, if relaxation is not complete at radii where par- 
ticles have near escape velocities, then the outer density 
profile slope may deviate from —4. 

4. DISCUSSION 

4.1. The approach to equilibrium 

In this paper, we used statistical mechanics to predict 
the energy distribution function of relaxed systems. A 
key feature of statistical mechanics approaches is that 
they deal only with the final equilibrium states, and de- 
rive these by equating them to the most probable or 
maximum entropy states. We stress that our theory is 
therefore limited to the description of the final state of 
self-gravitating collisionless systems. Because the final 




Figure 2. Statistical mechanical predictions for n (Equa- 
tion d 1 3 1) ^ - In the usual statistical mechanical approach, one identi- 
fies fi space with (x,v) phase space and n with / (upper panel). In 
the approach proposed in this paper, one identifies /i space with en- 
ergy space and n with N (lower panel). We have assumed </>o = 4. 
Different line styles signify different cutoffs: £ = (no cutoff, thin 
dashed line), f = exp(— 7) (thick solid cur ve, DARKexp), and dis- 
crete approach (thin solid line, Equation {TOjl ) . 

state is an equilibrium state, it is the result of full re- 
laxation and therefore does not rely on any additional 
assumptions. 

Implicit in our derivation is the assumption of equal 
a priori probabilities in state space. While this can 
be accomplished through efficient mixing (ergodicity), 
our theory does not address specific physical mecha- 
nism for attaining full m ixing. In a companion paper 
(| Williams fc Hiortbl I2010D . we use the Extended Sec- 
ondary Infall Model (ESIM) to test the DARKexp predic- 
tion, but stress that ESIM is a restricted physical model 
and is not equivalent to A-body simulations. 

A straightforward way to evaluate the applicability 
of physical mechanisms for relaxation is to compare 
their end states. For example, t he scattering model 
for vio lent relaxation introduced bv lSpergel fc HernquistJ 
(p92l) predicts N(E) oc [{E - $ )" 2 + C(E - 
$o) -1 ^ 2 ] -1 exp(— @E) which is clearly inconsistent with 
the DARKexp final state. On the other hand, the final 
state of ESIM halos is well fit with DARKexp. 

4.2. Collisional Versus Collisionless Systems 

Our prediction for N(e) applies to the end states of 
collisionless dynamics. Collisionlessness has several dif- 
ferent c onsequences, many of which can be used in one's 
theory. iLvnden-Belll (|1967f ) used the collisionless Boltz- 
mann equation to introduce an exclusion principle to ac- 
count for the incompressibility of the phase-space fluid. 
We use a different property; that collisionless dynamics 
implies the constancy of the energy per unit mass for 
each particle. Therefore, the collisionless nature of the 
problem is automatically included and there is no need 
for an exclusion principle. 
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In a system driven toward equilibrium by two-body 
relaxation, the situation is quite different. In this case, 
efficient relaxation implies that any particle can in prin- 
ciple end up anywhere in phase space and the usual par- 
tition of fi space is appropriate. In this case, takin g into 
accou nt low occupation numbers, one retrieves the iKingl 
(1966) f(E) which is rounded down at near escape en- 
ergies, while DARKexp's collisionless N(E) is rounded 
down at highly bound energies. 

The density profiles of the two are clearly different, 
but both are very good approximations to what one sees 
in globular clusters and simulated dark-matter haloes, 
respectively. 

4.3. Anisotropy 

In this paper we dealt exclusively with isotropic sys- 
tems for which the differential energy distribution is 
a function of energy only. The full simulated struc- 
tures, however, are definitely anisotropic, and exhibit 
a correlation between t he density slope and anisotropy 
(|Hansen fc Moordl2006T ). 

In principle, our statistical mechanics formalism can be 
extended to spherical non-isotropic systems, and we plan 
to do so in a future publication. In a fully developed the- 
ory, the distribution function must depend explicitly on 
angular momentum, es pecially in the outer parts (e.g., 
iStiavelli fc Bertirj I1987H . Whether the energy distribu- 
tion of the full theory will be significantly affected by 
the introduction of ang ular momentum is yet un clear; nu- 
merical experiments of lMacMillan et all (j2006L Figure 9) 
seem to indicate that the DARKexp form is still a good 
description of systems that underwent radial orbit insta- 
bility. 

For now we note that even though the DARKexp pre- 
diction may not be expected to be an excellent fit to 
simulations at all energies, it should apply to the cen- 
tral regions of dark-matter halos, as these are isotropic. 
Therefore, our explanation of the central density cusp of 
— 1 applies to simulated halos. 

The Dark Cosmology Centre is funded by the Danish 
National Research Foundation. L.L.R.W. is very grateful 
for the hospitality of the Dark Cosmology Centre at the 
University of Copenhagen and the Institute for Theoret- 
ical Physics at the University of Zurich, where she spent 
the Fall of 2009 and the Spring of 2010, respectively. 

APPENDIX 

LIMITING BEHAVIOR AT SMALL RADII FOR 
DIFFERENT EXPONENTIAL CUTOFFS 

One might accept that N(e) oc exp(— e) (Section 12. ip . 
but not the proposed form of the cutoff motivated by the 
low occupation numbers (Section 12. 2p . for example, be- 
cause our proposed transition from the integer to the con- 
tinuous form is more algebraic than physical. In this sec- 
tion we address possible alternatives to the cutoff shape. 

First, is a cutoff needed at all? If there is no cutoff, 
then the differential energy distribution is of the form 
Equation |8]). which would lead to a singularity in the 
central potential. To avoid that, one needs a cutoff. One 
possible form is to simply truncate N(s) at some finite 
central value; N(e) = C at s = 0o and N(e) — for 



e > 4>q. This would imply a — oo, and an infinite, pos- 
itive central density slope, i.e., a hole (Section [3- 1[) and 
divergent velocity dispersion, both of which are unphys- 
ical. 

Therefore, we conclude that a smooth cutoff is re- 
quired, i.e., that at some finite </>o, N(e = </>o) = 0. At 
energies near these highly bound energies, we can Taylor 
expand the non-truncated N(e). Using the abbreviation, 
x — 0o — exp(x) ~ 1 + x + x 2 /2 + • • • . Whatever the 
specific shape of the cutoff function, it too can be Tay- 
lor expanded, and then subtracted from that of exp(x). 
Three outcomes can result after this subtraction, depend- 
ing on the exponent of x of the surviving leading term, 

(1) x B , where B < 1 and x and higher terms cancel out; 

(2) x; and (3) x 2 , if the cutoff function goes exactly as 
1 + x. In case (1) a > 1, and the central density slope 
will be shallower than —1. Case (2) leads to our main 
prediction of the central density slope of —1, while case 

(3) is the only way to get density slopes steeper than 
— 1. In this case a = 0.5, and the central density slope is 
-1.5. 

Since mathematically it is possible to get central den- 
sity slopes shallower or steeper than —1, it is ultimately 
the physical arguments that will dictate the form of the 
cutoff in N(e) and hence the slope value. 
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